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Abstract 

Port-Hamiltonian systems result from port-based network modeling of physical systems and are an important example of 
passive state-space systems. In this paper, we develop the framework for model reduction of large-scale multi- input/multi- 
output port-Hamiltonian systems via tangential rational interpolation. The resulting reduced-order model not only is a rational 
tangential interpolant but also retains the port-Hamiltonian structure; hence is passive. This reduction methodology is described 

in both energy and co-encrgy system coordinates. We also introduce an H2-inspired algorithm for effectively choosing the 
interpolation points and tangential directions. The algorithm leads a reduced port-Hamiltonian model that satisfies a subset 
of H2-optimality conditions. We present several numerical examples that illustrate the effectiveness of the proposed method 
showing that it outperforms other existing techniques in both quality and numerical efficiency. 
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1 Introduction 

Port-based network modeling [37] of physical systems ex- 
ploits the common circumstance that the system under 
study is decomposable into (possibly many) subsystems 
which are interconnected through (vector) pairs of vari- 
ables, whose product gives the power exchanged among 
subsystems. This approach is especially useful for multi- 
physics systems, where the subsystems may describe dif- 
ferent categories of physical phenomena (e.g, mechan- 
ical, electrical, or hydraulic). This approach leads di- 
rectly to port-Hamiltonian system representations (see 
[38,40,39,37] and references therein) encoding structural 
properties related to the manner in which energy is dis- 
tributed and flows through the system. 
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Quite apart from underlying port-Hamiltonian struc- 
ture, models of complex physical systems often involve 
discretized systems of coupled partial differential equa- 
tions which lead, in turn, to dynamic models having 
very large state-space dimension. This creates a need for 
model reduction methods capable of producing (com- 
paratively) low dimension surrogate models that are 
able to mimic closely the original system's input / output 
map. When the system of interest additionally has port- 
Hamiltonian structure, then it is highly desirable for 
the reduced model to preserve port-Hamiltonian struc- 
ture so as to maintain a variety of system properties, 
such as energy conservation, passivity, and existence of 
conservation laws. 

Preservation of port-Hamiltonian structure in the reduc- 
tion of large scale port-Hamiltonian systems has been 
considered in [30,28,32] using Gramian-based methods; 
we briefly review one such approach in Section 3.1. For 
the special case of a single- input /single-output (SISO) 
system, the use of (rational) Krylov methods was em- 
ployed in [26,15,27,44]. In particular, [26,27,44] all deal 
with system interpolation at single points only. In this 
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work, wc develop multi-point rational tangential in- 
terpolation of multi-input/multi-output (MIMO) port- 
Hamiltonian systems. Preservation of port-Hamiltonian 
structure by reducing the underlying full-order Dirac 
structure was presented in [31,41]. A perturbation ap- 
proach is considered in [19,18]. Sec [25] for an overview 
of recent port-Hamiltonian model reduction methods. 

For general MIMO dynamical systems, interpolatory 
model reduction methods produce reduced-order mod- 
els whose transfer functions interpolate the original sys- 
tem transfer function at selected points in the complex 
(frequency) plane along selected input and output di- 
rections. The main implementation cost involves solving 
(typically sparse) linear systems, which gives a signif- 
icant advantage in large-scale settings over competing 
Gramian-based methods (such as balanced truncation) 
that must contend with a variety of large-scale dense 
matrix operations. Indeed, a Schur decomposition is re- 
quired for computing an exact Gramian although recent 
advances using ADI-type algorithms have allowed the it- 
erative computation of approximate Gramians for large- 
scale systems; see, for example, [24,6,16,36,34,20,9,33] 
and references therein. 

Until recently, there were no systematic strategies for 
selecting interpolation points and directions, but this 
has largely been resolved by Gugercin et al. [12,13,14] 
where an interpolation point /tangent direction selection 
strategy leading to 'H2-optimal reduced order models has 
been proposed. See [42,7,22] for related work and [2] for 
a recent survey. 

The goal of this work is to demonstrate that interpola- 
tory model reduction techniques for linear state-space 
systems can be applied to MIMO port-Hamiltonian sys- 
tems in such a way as to preserve the port-Hamiltonian 
structure in the reduced models, and preserve, as a con- 
sequence, passivity as well. Moreover, we introduce a 
numerically efficient 'H2-based algorithm for structure- 
preserving model reduction of port-Hamiltonian systems 
that produces high quality reduced order models in the 
general MIMO case. Numerical examples are presented 
to illustrate the effectiveness of the proposed method. 

The rest of the paper is organized as follows: In Section 2, 
wc review the solution to the rational tangential interpo- 
lation problem for general linear MIMO systems. Brief 
theory on port-Hamiltonian systems is given in Section 
3. Structure-preserving interpolatory model reduction 
of port-Hamiltonian systems in different coordinates is 
considered in Section 4 where we show theoretical equiv- 
alence of interpolatory reduction methods using different 
coordinate representations of port-Hamiltonian systems 
and discuss which is most robust and numerically ef- 
fective. 'H2-based model reduction for port-Hamiltonian 
systems together with the proposed algorithm is pre- 
sented in Section 5 followed by numerical examples in 
Section 6. 



2 Interpolatory Model Reduction 

Let G(s) be a dynamical system with a state-space re- 
alization given as 



G(a): 



Ex(i) = Ax{t) + Bu{t), 
y(i) = Cx(t), 



(1) 



where A,E e M"""", E is nonsingular, B e M"""™, 
and C G RP''". In (1), for each t we have: x(t) e M", 
u{t) e M™ and y{t) e denoting, respectively, the 
state, the input, and the output of G(s). By taking the 
Laplace transform of (1), we obtain the associated trans- 
fer function G(s) = C(sE — A)~-^B. Following usual 
conventions, the underlying dynamical system and its 
transfer function will both be denoted by G(s). 

We wish to produce a surrogate dynamical system Gi.(s) 
of much smaller order with a similar state-space form 



Oris) 



Er±r{t) = ArXr(t) -|- BrU{t) 
Yrit) = C,X,,(t), 



(2) 



where A^, E^ e W'''', e W""", and € W'' with 
r <^n such that the reduced-order system output y,. (t) 
approximates the original system output, y{t), with high 
fidelity, as measured in an appropriate sense. The ap- 
propriate sense here will be with respect to the H2 norm 
and will be discussed in Sections 2.2 and 5. 

We construct reduced order models Gr(s) via Petrov- 
Galerkin projections. We choose two r-dimensional sub- 
spaces Vr and Wr and associated basis matrices G 

M"'''' and e K"'"' (such that Vr = Range(Vr) and 
yVr = Range(Wi.)). System dynamics arc approximated 
as follows: we approximate the full-order state x[t) as 
x(i) « VrXr(t) and force a Petrov-Galerkin condition as 

(EVriCr - AVrXr - B u) = 0, = CYrXr, (3) 

leading to a reduced-order model as in (2) with 



E^ = W^EV^, B^ = W^B, 
r = WTAV., and = CV.. 



(4) 



The quality of the reduced model depends solely on the 
selection of the two subspaces Vr and Wr (or equivalently 
on the choices of and W^.). We choose Vr and 
to force interpolation. 

2. 1 Interpolatory projections 

The construction of reduced order models with inter- 
polatory projections was initially introduced by Skelton 
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et. al. in [8,46,47] and later put into a robust numeri- 
cal framework by Grimme [11]. This problem framework 
was recently adapted to MIMO dynamical systems of 
the form (1) by Gallivan et al. [10]. 

Beattie and Gugercin in [4] later extended this frame- 
work to a much larger class of transfer functions, namely 
those having a generalized coprime factorization of the 
form H(s) = C(s)K(.s)-iB(s) with B(s), C(s), and 
K(s) given as meromorphic matrix-valued functions. 
Given a full-order system (1) and interpolation points 
{skYk=i ^ wifh corresponding (right) tangent di- 
rections {bfcjJJ^j e C™, the tangential interpolation 
problem seeks a reduced-order system G,.(,s) that inter- 
polates G(s) at the selected points along the selected 
directions; i.e.. 



G{si)bi = Gr(sj)bi, 



for i — \, - ■ ■ 



(5) 



Analogous left tangential interpolation conditions may 
be considered, however (5) will suffice for our purposes. 
Conditions forcing (5) to be satisfied by a reduced system 
of the form (4) are provided by the following theorem 
(see [10]). 

Theorem 1 Suppose G{s) ~ C(sE — A)^-^B. Given 
a set of distinct interpolation points {skYk=i ^'^'^ right 
tangent directions {bfe}^^;^, define Yr € C"^*" as 



provided that sE — A and sE^ — A,, are invertible where 
G^^\s) denotes the£^^ derivative o/G(s). By combining 
this result with Theorem 1, one can match different inter- 
polation orders for different interpolation points ( analo- 
gous to generalized Hermite interpolation). For details; 
see, for example, [10,2]. 

Remark 4 Notice that in Theorem 1 and Remark 3, 

what guarantees the interpolation is the range of the ma- 
trix Vj., not the specific choice of basis in ( 6). Hence, one 
may substitute with any matrix having the same 
range. In other words, for any satisfying V,. = V^L 
where L G W^"^ invertible, the interpolation property 
still holds true. This is a, simple consequence of the fact 
that the basis change L from to amounts to a 
state-space transformation in the reduced model. We will 
make use of this property in discussing the effect of dif- 
ferent coordinate representations for port-Hamiltonian 
systems. When the interpolation points occur in com- 
plex conjugate pairs ( as always occurs in practice ), then 
the columns o/ also occur in complex conjugate pairs 
and may be replaced with a real basis, instead. We write 
Vr = |vi, . . . , Wr\ to represent that a real basis for 
is chosen. Notice then that will also be a real matrix 
which then leads to real reduced order quantities as well. 

2.2 Interpolatory optimal approximation 

The selection of r interpolation points {si}^^^ £ C and 



Vr = [(siE — A) ^Bbi, ••• , (srE — A) ^Bb^] (6) corresponding tangent directions {b^} 



Then for any G 

Gr{s) = Cr(sEr - A^ 

(5) provided that SjE - 
fori = l,...,r. 



n"^*^, the reduced order system 
)~^'Br defined as in (4) satisfies 
A and s,;Er — A^ are invertible 



Remark 2 Attempting interpolation of the full transfer 
function matrix (in all input and output directions) gen- 
erally inflates the reduced system order by a factor of 
m (the input dimension). However, this is neither nec- 
essary nor desirable. Effective reduced order models, in- 
deed, 'H2 -optimal reduced order models, can be generated 
by forcing interpolation along selected tangent directions. 

Remark 3 The result of Theorem 1 generalizes to 
higher-order interpolation (analogous to generalized 
Hermite interpolation) as follows: For a point s € C and 

tangent direction b, suppose 



((sE - A)-^ E)*" ^ (sE - A)-^Bb € Range(V,.) 
fork = l,...,N. 



e C" com- 
pletely determines projecting subspaces and so is the 
sole factor determining the quality of a reduced-order 
models. Until recently only heuristics were available to 
guide this process and the lack of a systematic selection 
process for interpolation points and tangent directions 
had been a key drawback to interpolatory model reduc- 
tion. However, Gugercin et al. [14] introduced a frame- 
work for determining optimal interpolation points and 
tangent directions to find optimal 'H2 reduced-order ap- 
proximations, and also proposed an algorithm, the It- 
erative Rational Krylov Algorithm (IRKA) to compute 
the associated reduced-order models. In this section we 
briefly review the interpolatory optimal H2 problem for 
the general linear systems G(s) = C(sE — A)~^B with- 
out structure. 

Let Ai{r) denote the set of reduced order models having 
state-space dimension r as in (2). Given a full order sys- 
tem G(s) = C(sE - A)^^B, the optimal-'H2 reduced- 
order approximation to G(s) of order r is a solution to 



(7) 



G — Gr 



Then, for any e C"^'', the reduced order system 
Gr{s) defined as in (4) satisfies 

GW (s)b = GW {s)b, for £ = 0,--- ,N-1, (8) 



where 



iiGik ■.= [^j_jGi^)r^d^ 



■H2 



1/2 



(9) 
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There are a variety of methods that have been de- 
veloped to solve (9). They can be categorized as 
Lyapunov-based methods such as [45,35,17,21,43,48] 
and interpolation-based optimal methods such as 
[23,14,13,42,7,12,22,3,5]. Although both frameworks are 
theoretically equivalent [14], interpolation-based opti- 
mal 'H2 methods carry important computational ad- 
vantages. Our focus here will be on interpolation-based 
approaches. 

The optimization problem (9) is nonconvex, so obtain- 
ing a global minimizer is difficult in the best of circum- 
stances. Typically local minimizers are sought instead 
and as a practical matter, the usual approach is to find 
reduced-order models satisfying the first-order necessary 
conditions of ■H2-optimality. 

Interpolation-based ?^2-optimality conditions for SISO 
systems were introduced by Meier and Luenberger [23]. 
Interpolatory 'H2-optimality conditions for MIMO sys- 
tems have recently been developed by [14,7,42] . 

Theorem 5 Given a full-order system G(s) — C(sE — 
A)-iB, letGr{s) = -^Zx '^-^T bestr^"^ order 

approximation of G with respect to the 'H2 norm. Then 



(a) G(-Afc)bfc = Gr(-Afc)bfc, 

(h) c^G(-Afe) = Cfc Gr(-Afc), and 
(c) clG'{-Xk)hk = 4G;(-Afc)bfc 



(10) 

(11) 
(12) 



fork = 1, 2, r. 



Theorem 5 asserts that any solution to the optimal- 7^2 
problem (9) must be a bitangential Hermite intcrpolant 
to G(s) at the mirror images of the reduced system poles. 
This woiild be a straightforward construction, if the re- 
duced system poles and residues were known a priori. 
However, they arc not and the Iterative Rational Krylov 
Algorithm, IRKA [14] resolves the problem by iteratively 
correcting the interpolation points by refiecting the cur- 
rent reduced-order system poles across the imaginary 
axis to become the next set of interpolation points and 
the directions. The next tangent directions are residue 
directions taken from the current reduced model. Upon 
convergence, the necessary conditions of Theorem 5 are 
met. IRKA corresponds to applying interpolatory model 
reduction iteratively using 



Wr 



[(siE - A)-iBbi, • • • , {srB - A)-iBb,] 
= [(siE - A)-'^C^ci, • • • , {srE - A)-^C^c.] 



(13) 



where Sj, bj and Cj are, respectively, the interpolation 
points, and right and left tangent directions at step k of 
IRKA. Note that the Hermite bitangential interpolation 
conditions (10)-(12) enforce a choice on as in (13) 
in contrast to Theorem 1 where could be chosen ar- 
bitrarily since bitangential interpolation is not required. 



We revisit this issue in Section 5. For details on IRKA, 
see [14]. 



3 Linear port-Hamiltonian Systems 



In the absence of algebraic constraints, linear port- 
Hamiltonian systems take the following form ([38,37,28]) 



G(s): 



X = (J-R)Qx-FBu, 
y = B^Qx, 



(14) 



where J = -J^, Q = Q"^, and R = R'^ > 0. Q is 

the energy matrix; R is the dissipation matrix; J and B 
determine the interconnection structure of the system. 
H{x) = ix-^Qx defines the total energy (the Hamilto- 
nian). For all cases of interest, -ff(x) > 0. Note that the 
system (14) has the form (1) with E = I, A = (J-R)Q, 
and C = B^Q. 



The state variables x G R" are called energy variables, 
since the total energy i?(x) is expressed as a function 
of these variables, u, y G M™ are called power variables, 
since the scalar product u-^y equals the power supplied 
to the system. Note that 



dt 



Qx = u^y-x^QRQx < u^y 



— that is, port-Hamiltonian systems are passive and the 
change in total system energy is bounded by the work 
done on the system. 



We will assume throughout that Q is positive definite. In 

the case of zero dissipation (R = 0), the full-order model 
G(.s) of (14) is both lossless and stable. For positive 
semidefinite R, G(,s) could be either stable or asymptot- 
ically stable in general; and for R strictly positive defi- 
nite, the system is asymptotically stable. We will assume 
that G(s) is merely stable except in Section 5 where we 
assume asymptotic stability in order to assure a bounded 
•^2 norm. 



Example 1 Consider the Ladder Network illustrated in 
Fig. 1, withCi, C'2, Li, L2, Ri, R2, R3 being, respectively, 
the capacitances, inductances, and resistances of ideal- 
ized linear circuit elements described in the figure. The 
port-Hamiltonian representation of this physical system 
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has the form (14) with 



Ri Li,(/>i 



i2, 02 



U2 = U~ 



-1 -1 1 

J = tridiagi 
11-1 



R = diag(0,i?i,0,i?2 + i?3), and B 



The A matrix is 



A = 



1 


1 



(15) 













1 ^-1 





LI' 










-c,-^ - 


(i?2 + Ri)L2 ^ 


space 


vector X is 


gitien as 



gi 01 92 02 



wzi/i gi, 92 being the charges on the capacitors Ci, C2 and 
01,02 feemj </ie fluxes over the inductors Li,L2 corre- 
spondingly. The inputs of the system, {ui,M2} are given 
by the current I on the left side and the voltage U on the 
right side of the Ladder Network. The port-Hamiltonian 
outputs {^1,1/2} are the voltage over the first capacitor 
Uci and the current through the second inductor Ij^^ . 

The system matrices A, B, C, andE of (1) follow directly 
from writing the linear input-state differential equation 
for this system. Q may be derived from the Hamiltonian 
-ff (x) = ix^Qx. Once A and Q are known, it is easy 
to derive the dissipation matrix, R, and the structure 
matrix, J, corresponding to the Kirchhoff laws, such that 
A = (J — R)Q. The port-Hamiltonian output matrix C 
is given as 



C = B^Q = 



^000 

000' 



L2 J 



(16) 



Fig. 1. Ladder Network 

The coordinate transformation [37,25] between energy 
coordinates, x, and co-energy coordinates, e, is given by 



Qx. 



(18) 



Example 2 (continued). The co-energy state vector, e, 
for the Ladder Network from Example 1 is given as 



Uci Ili 



with , being the voltages on the capacitors Ci , C2 
and , being the currents through the inductors 
Li,L2, respectively. 

3. 1 Balancing for port-Hamiltonian systems 

To make the presentation self-contained, we review a re- 
cent structure-preserving, balancing-based model reduc- 
tion method for port-Hamiltonian systems, the effort- 
constraint method ([31,41,25,32,29]), which will be used 
for comparisons in our numerical examples. 

Consider a full order port-Hamiltonian system realized 
as in (14) with respect to energy coordinates. Consider 
the associated balancing transformation, T;,, defined in 
the usual way (see [1] for a complete treatment): 
simultaneously diagonalizes the observability Gramian, 
Go, and the controllability Gramian, Gc, so that 



3 -^h 



diag((Ti, 0-2, ■ • • , cr„). 



where cti, (T2, • ■ • , c„ are the Hankel singular values of 
the system G(s). This balancing transformation is the 
same that is used in the context of regular balanced trun- 
cation. 

The state-space transformation associated with balanc- 
ing (as with any other linear coordinate transforma- 
tion), will preserve port-Hamiltonian structure of the 
system (14). Indeed, if we define balanced coordinates, as 
Xf, = TfoX, then we have 



Recall from [37,25,28] that the port-Hamiltonian system 
(14) may be represented in co- energy coordinates as 



e = Q(J R)e + QBu, 



(17) 



Xfe = (Jb - R;,)QbXb + Bf,u, 
y = (Bb)^QbXb, 



(19) 



where T^RT^ = Rf, = R^ > is the dissipation ma- 
trix, TbJT?" — 3h — —JT is the structure matrix and 



5 



'^QT, ^ = Qb = > is the energy matrix in bal- 



^-^b — — ^b 

anced coordinates. In this case, B(, 



ThB. 



4.1 Interpolatory projection with respect to energy co- 
ordinates 



Now split the state-space vector, X{,, into dominant and 

Xfcl 



codominant components: Xb = 



X62 



where Xf,i e 



Regular balanced truncation proceeds by truncating the 
codominant state-space components, in effect forcing the 
system to evolve under the constraint X(,2 = 0. This 
will destroy port-Hamiltonian structure in the remaining 
reduced order system that determines the evolution of 

Xbl. 

In contrast to this, the effort- constraint method produces 
the following reduced order port-Hamiltonian system 



X6I = (J6II — R-61l)S6Xbi -I- BfciU, 
Yr = (B6i)^S6Xfti, 



(20) 



where S;, = Qbii - Qbi2(Qfc22)~^Qf)2i is the Schur com- 
plement of the energy matrix Qf, in balanced coordi- 
nates; the other matrices are of corresponding dimen- 
sion. The relationship of the effort-constraint method to 
the reduction of the underlying (full-order) Dirac struc- 
ture is discussed in [31,41,25]. A more direct approach 
is given in [29] ; another approach using scattering coor- 
dinates can be found in [32] . 

Remark 6 The effort- constraint method can he formu- 
lated with Petrov-Galerkin projections [31,41,25]. Even 

though a balancing transformation is used in develop- 
ing the effort- constraint method as expressed in (20), the 
method is not equivalent to the usual balanced truncation 
method. Note that balanced truncation does not preserve 
port-Hamiltonian structure. For more details, see [25]. 



We first show how to achieve interpolation and preserva- 
tion of port-Hamiltonian structure simultaneously in the 
original energy coordinate representation. The choice of 
Wr plays a crucial role. 

Theorem 7 Suppose G(s) is a linear port-Hamiltonian 
system, as in (14). Let {si]\^i C C be a set ofr distinct 
interpolation points with corresponding tangent direc- 
tions {bi}l^i C C", both sets being closed under con- 
jugation. Construct as in (6) using A = (J — R)Q 
and E = I, so that 

Vr = [(siI-(J-R)Q)-iBbi, . . . , (s.I-(J-R)Q)-iBb.]. 

LetM. e C^'' 6e any nonsingular matrix such thatVr 
\ \ I is real. 

Define = QV^(V^QVr)"^ and then 



= WjRW. 



and Br = WTB. 



(21) 



Then, the reduced-order model 



Gr{s) : 



(22) 



is port-Hamiltonian, passive, and 



G{si)bi = Gr{si)bi 



fori = l,--- ,r. 



That is, Gr(s) interpolates G(s) at {si}l_i along the 
tangent directions {bi}[^j. 



4 Interpolatory model 
Hamiltonian systems 



reduction of port- Proof: In Ught of Remark 4, there is at least one (and 

hence, an infinite number) of possible M e satis- 



Now we consider how to create an interpolatory reduced- 
order model that retains port-Hamiltonian structure in- 
herited from the original system, thus guaranteeing that 
it is both stable and passive. 

Note that Theorem 1 may be applied to a port- 
Hamiltonian system with an arbitrary choice of W^.. 
As a result, will have the form W^(J — R)QVr, 
which can be no longer represented as (J^ — Rr)Qr 
with a skew-symmetric J^, and symmetric positive 
semi-definite Rr and positive definite matrices. Be- 
low, we develop several approaches that will choose 
carefully to allow the structure preservation. 



tying the conditions of the Theorem and so, V,. e M"^''. 
Trivially, one may observe that the reduced-order sys- 
tem (22) is real and retains port-Hamiltonian structure 
with a positive definite Q^, and so must be passive. To 
show that Gr(s) interpolates G(s) note, as before, that 
the port-Hamiltonian system (14) has the standard form 
(1) with E = I, A = (J-R)Q, and C = B^Q. An inter- 
polatory reduced model may be defined as in (2) using 

Yr = Yr andJW^ = W^. We then have E^ = W^V^ = 
Ir and A^ = W^AV^ = W^(J - R)QV^. From Theo- 
rem 1, this model interpolates G(s) at {si}^^^ along the 
tangent directions {bi}[_;^ as required but it is not obvi- 
ous that this is the same reduced system Gr{s) that we 
have defined above. 
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Note that W^V^ is a (skew) projection onto Range(QVr) 
and so QV^ = WrV^QV^. Thus, we have 

=W^(J - R)QV^ 
=W^(J - R)WrV^QV^ 

and 



=CV^ = B^QV^ 

which verifies that the port-Hamiltonian reduced sys- 
tem, Gr{s), interpolates G(s) at {si}^^^ along {bj}^^^ 
as required. □. 

Remark 8 Theorem 7 can be generalized to include 

higher-order (generalized, Hermite) interpolation fol- 
lowing ideas described in Remark 3. Indeed, given an 
interpolation point s and tangent direction b, augment 
Vr so that (7) holds. The remaining part of the con- 
struction of Theorem 7 proceeds unchanged and the re- 
sulting reduced- order model will retain port-Hamiltonian 
structure (hence be both stable and passive) and will 
interpolate higher-order derivatives o/G(s) as in (8). 



4-2 The influence of state-space transformations 

Let T S M"^" be an arbitrary invertible matrix rep- 
resenting a state-space transformation x = Tx. As 
one may recall from the discussion of Section 3.1, 
port-Hamiltonian structure will be maintained under 
state-space transformations: 

X = (J - R)Qx-|-Bu, X = (J- R)Qx-|-Bu, 

y = B^Qx. y = B^Qx. 

with transformed quantities 

J = TJT^, R = TRT^, 
Q = T-^QT-i, and B = TB. 



Proof: Define primitive interpolatory basis matrices for 
each realization as 

V. = [(siI-(J-R)Q)-iBbi,...,(s.I-(J-R)Q)-iBb.] 
and 

V. = [(siI-(J-R)Q)-iBbi,...,(s.I-£J-R)Q)-iBb.l. 

Elementary manipulations show that = TV^.. Note 
that any M € C'^^'^ for which V,. = V^M is real makes 

Yr = VrM real as well (and vice versa). Then we define, 
following the construction of Theorem 7, = QV^ 

and Wr = QVr Observe that = T-'^W^ and 
= TV^. So, 



= V^QV^ 
R, = WjRW, 



v; QV^ = 



wJrw, 



and B^ = W?^B = W^'b = B^. 



□ 



Remark 10 One may conclude from Theorem 9 that 
prior to calculating a reduced- order interpolatory model, 
there may be little advantage in first applying a state- 
space transformation (e.g., a balancing transformation, 
Tb, as described in Section 3.1 or transforming to co- 
energy coordinates with T Q as m (18)). Indeed, 
choosing a realization that exhibits advantageous spar- 
sity patterns facilitating the linear system solves neces- 
sary to produce V,. 2ijill likely be the most effective choice. 
If we choose T to satisfy T"^T = Q (for example, if 
T is a Cholesky factor for Q ) then, with respect to the 
new yi- coordinates (called "scaled energy coordinates"), 

we find that Q = I and so — V^. Notice that in this 
case, scaled energy coordinates and scaled co- energy co- 
ordinates become identical. Notwithstanding these sim- 
plifications, unless the original Q is diagonal (so that the 
transformation to scaled energy coordinates is cheap), 
there appears to be little justification for global coordinate 
transformations preceding the construction of an inter- 
polatory reduced order model. 

5 H2-based reduction of port-Hamiltonian sys- 
tems 



Theorem 9 Suppose G(s) is a port-Hamiltonian sys- 
tem as in (14), with two different port-Hamiltonian 
realizations connected via a state-space transformation, 
T, as above. Given a set of r distinct interpolation 
points {si}l^i G C with corresponding tangent directions 
{^i}i=i G (closed under conjugation), then interpo- 
latory reduced-order port-Hamiltonian models produced 
from either realization via Theorem 7 will be identical to 
one another. 



Inspired by the optimal '^2 model reduction method de- 
scribed in Section 2.2, we propose here an algorithm that 
produces high quality reduced-order port-Hamiltonian 
models of port-Hamiltonian systems. For the 'H2 norm 
to be well defined, we assume that the full-order model 
G(s) is asymptotically stable. 

Even though the optimality conditions (10)-(12) can be 
satisfied for generic model reduction settings starting 
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with a general system G(s) — C(sE — A)^^B and re- 
ducing to Gr{s) = Cr(sEr — Ar)~^Br without any re- 
striction on structure, it will generally not be possible 
to satisfy all optimality conditions while also retaining 
port-Hamiltonian structure; there are not enough de- 
grees of freedom to achieve both goals simultaneously. 
Note that (10)-(12) corresponds to r{pm, + 1) condi- 
tions — precisely the number of degrees of freedom in 

Gr{s) = Cris'Er ~ A^)-lB^ = ^[^^ ^"^i^f . This 

representation is a general partial fraction expansion 
and need not correspond to a port-Hamiltonian system. 
Adding port-Hamiltonian structure on top of (10)-(12) 
creates an overdetermined set of optimality conditions 
that cannot typically be satisfied. Alternatively, note 
that the optimality conditions (10)-(12) requires choos- 
ing Wr in a specific way as in (13) unlike the way it was 
chosen in Theorem 7 in order to preserve structure. See 
also Remark 13. 

We chose to enforce only a subset of first-order optimal- 
ity conditions and use the remaining degrees of freedom 
to preserve structure. We enforce (10) with a modifed 
version of IRK A that maintains port-Hamiltonian struc- 
ture at every step. A description of our proposed algo- 
rithm is given as Algorithm 1. 

Algorithm 1 IRKA for MIMO port-Hamiltonian 
systems (IRKA-PH): 

(1) LetG{s) = B^Q(sI-(J-R)Q)-iB be as in (14). 

(2) Choose initial interpolation points {si, . . . , Sr} and 
tangential direetions {bi, . . . , br}; both sets closed 
under conjugation. 

(3) Construct 

V, = [(siI-(J-R)Q)-iBbi,..., 

...,(sj- (J-R)Q)-iBb^]. 

(4) Choose a nonsingular matrix M e C'"^'" such that 
Vr = VrM is real. 

(5) CalculateWr = Q'Vr(VjCiVr)~'^ 

(6) until convergence 

(a) 3r = W^JW^, Rr = W^jaW^, 

= V^QV^ and Br = WjB. 

(b) A,. = (J^ - R^)Q,. 

(c) Compute ArXi = XiXi, y*A^ = A^y* 

with y*Xj = 6ij for left and right eigenvectors 

y* and associated with Xi . 

(d) Si < Xi and bf ^ y^ for i = 1, . . . , r. 

(e) Compute = |(.sil - (J - R)Q)-iBbi, • ■ • , 

...,(sJ-(J-R)Q)-iBb^l. 

(f) Choose a nonsingular matrix M G g^ch 
that Vr = V^M is real. 

(g) Ca/cwtoeW^ = QV^(V^QV^)-i 

(7) The final reduced model is given by 

Jr = W^JW^, R^ = W^RW^, B^ = W^B, 
Qr = V'^QVr, and Cr = B^Q^. 



Theorem 11 LetG{s) = B^Q(sI- (J-R)Q)-iB bea 
port-Hamiltonian system as in (14). Suppose IRKA-PH 
as described in Algorithm 1 converges to a reduced model 

Gr(s) = ^^ibf. (23) 
i=l Xi 

Then Gr{s) is port-Hamiltonian, asymptotically stable, 
and passive. Moreover, G,.(s) satisfies the necessary con- 
dition (10), for 1-12 optimality: G,.(,s) interpolates G(s) 
at —Xi along the tangent directions bj for i = 1, . . . ,r. 

Proof: The port-Hamiltonian structure, and conse- 
quently passivity, are direct consequences of the con- 
struction of Gr.(s) as shown in Theorem 7. The 7^2- 

based tangential interpolation property results from the 
way the interpolation points Si and the tangential direc- 
tions bj are corrected in Steps 6-c) and 6-d) throughout 
the iteration. Hence upon convergence; Si = and bi 
is obtained from the residue of Gr(s) corresponding to 
Aj as desired. 

To prove asymptotic stability we proceed as follows: 
Due to Theorem 9, withoiit loss of generality we can 
assume that Q = I. Choose M in IRKA-PH so that 
Yr is orthogonal, i.e. V^V^ = Ir Hence, we obtain 
Wr = Vr-. Then, the reduced-order quantities are given 
by the one-sided projection A^ = VjAV^, and B^ = 
= V^B. Let Ar = X^AX"^ be the eigenvalue de- 
composition of Ar upon convergence oi IRKA-PH where 
A = diag(Ai, . . . , Ai). Note that upon convergence, the 
tangential directions bi are obtained by transposing the 
rows of the reduccd-B matrix of Gr(s) represented in 
the modal form. Define = [bi, . . . , b^]^. Hence, upon 
completion of IRKA-PH, we have 

A^ = V^AV^ = X^AX^i and B^ = V^B = X^F^ 

From the construction of Vj. in Step (6)-(e) of IRKA- 
PH and from the fact that Si ~ —Xi upon convergence, 
the i^^ column of V,- satisfies (— Ajl — A)vj = Bbj. Con- 
sequently, Vr solves a special Sylvester equation 

AV^ + V^A^ + BF^ = (24) 

where we used the fact A ~ A^. Multiply (24) by X^ 
from right and use = V^M to obtain 

AV^K^ V^K^A^ + BBj = (25) 

where = M~^X^ e C^'' is non-singular. Then mul- 
tiplying (25) by V"^ from left yields 

A^K^ -\- K^A^ -\- B^B^ = (26) 
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Since Gr(s) is passive, we know that A,, has no eigenval- 
ues in the open right- half plane. Hence, to prove asymp- 
totic stability of Gr(s), we need to further show that 
has no eigenvalues on the imaginary axis either. Assume 
to the contrary that it has, i.e. 

z*Ar = z*A, with A = jw. 

Note that, A^z^ = Az^ = -Az^. Multiply (26) by 2.* 
from left and by z^ from right to obtain 

B^z^ = 0. (27) 

Then, multiplying (25), from right, by z^ and using (27) 
leads to 

AVrKrZ-r = V^K^ZrA. 

However, this means that V^K^z^ is an eigenvector of 
A with the corresponding eigenvalue A = jw, contra- 
dicting asymptotically stability of G(s). Therefore, A^ 
cannot have a purely imaginary eigenvalue and Gr{s) 
is asymptotically stable. □ 
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Fig. 2. Mass-Spring-Dampcr system 
optimality, (11) and (12), will be satisfied if in addition: 

Range[(AiI + (J - R)Q)-^Bbi, . . . , 

(A4 + (J-R)Q)-iBb.] = 

Range[(AiI + (J - R)^Q)-iBci, • ■ • , 

(A4 + (J-R)'^Q)-iBcr], 

where Aj are reduced system poles and bi, Cj are the cor- 
responding (vector) residues in the expansion (23) for 
i = 1, . . . ,r. The condition (29) allows retention ofport- 
Hamiltonian structure while enforcing the'H2- optimal hi- 
tangential interpolation conditions (10)-(12). Typically 
(29) will not hold, however the special case when J = 
will lead to (29) with hi — Ci implying satisfaction of 14.2 
optimality conditions for the reduced system. 



The convergence behavior of IRK A has been studied in 
detail in [14] ; see, also [2] . In the vast majority of cases 
IRKA converges rapidly to the desired first-order con- 
ditions independent of initialization. Effective initializa- 
tion strategies have been proposed in [14] as well, al- 
though random initialization often performs very well. 
We illustrate the robustness of IRKA-PH with regard 
to initialization in Section 6. 

Remark 12 Let'PH{r) denote the set of port- Hamiltonian 
systems with state-space dimension r as in (22) and 
consider the problem: 



IG-G, 



\H2 



mm 



G- G. 



■H2 



(28) 



where G{s) andGr{s) are both port- Hamiltonian. This is 
the problem, that we would prefer to solve and it is a topic 
of current research. Note that although IRKA-PH gen- 
erates a reduced-order port-Hamiltonian model, Gr{s), 
satisfying a first-order necessary condition for 'H2 opti- 
mality, this condition is not a necessary condition for 
Gj,{s) to solve (28), that is, for Gr{s) to be the opti- 
mal reduced-order port-Hamiltonian system approxima- 
tion to the port-Hamiltonian system G(s). Nonetheless, 
we find that reduced-order models produced by IRKA-PH 
have much superior 7^2 performance compared to other 
approaches. This is illustrated in Section 6. 

Remark 13 IRKA-PH, as described in Algorithm 1, 
produces a reduced order port-Hamiltonian model Gr{s), 
which satisfies the first-order necessary condition of 712- 
optimality (10). The remaining two conditions for 712- 



6 Numerical Examples 

In this section, we illustrate the theoretical discussion on 
three port-Hamiltonian systems. The first two systems 
are of modest dimension so that we can compute all the 
norms (such as H,2 and Hoo error norms) explicitly for 
several r values for many different methods to compare 
them to the fullest extent. Then, to illustrate that the 
proposed method can be easily applied in the large-scale 
settings as intended, we use a large-scale model as the 
third example. 

6. 1 MIMO Mass-Spring-Damper system 

The full-order model is the the Mass-Spring-Damper sys- 
tem shown in Fig. 2 with masses mi, spring constants kt 
and damping constants q > 0, i = 1, . . . , n/2. qi is the 
displacement of the mass m^. The inputs ui,U2 are the 
external forces applied to the first two masses mx,m2. 
The port-Hamiltonian outputs yi , 2/2 are the velocities 
of masses mi,TO2- The state variables are as follows: xi 
is the displacement qi of the first mass mi , X2 is the mo- 
mentum pi of the first mass mi , X3 is the displacement 
q2 of the second mass m2, X4, is the momentum p2 of the 
second mass m2, etc. 

A minimal realization of this port-Hamiltonian system 
for order n = 6 corresponding to three masses, three 
springs and three dampers is 
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Then the A matrix can be obtained as 



A = (J - R)Q = 






1 

mi 














-ki 


_ 

mi 


ki 




















1 

m2 








ki 





—ki - k2 


_ _£2_ 

m2 


k2 




















1 

ma 








k2 





-k2 - ks 


C3 

ms - 



Adding another mass with a spring and a damper would 
increase the dimension of the system by two. This wiU 
lead to a zero entry in the (n — l,n — 1) position and 
an entry of — c„/2/™n/2 in t^^c (n,n) position. The su- 
perdiagonal of A will have kn/2-1 in the (n — 2, n — 1) 
position and l/m„/2 in the (n — 1, n) position. The sub- 
diagonal of A will have in the (n — 1, n — 2) position 
and —kn/2-i — kn/2 in the (n, n — 1) position. Addition- 
ally A will have kn/2-1 in the (n, n — 3) position. 

We used a 100-dimensional Mass-Spring-Damper sys- 
tem with TUi — 4, ki — 4, and q = 1. We compare three 
different methods: (1) The proposed method in Algo- 
rithm 1 denoted by IRKA-PH (2) The effort-constraint 
method (20) of ([31,41,25,32,29]) denoted by Eff. Bal, 
and (3) One step interpolatory model reduction without 
the 7^2 iteration, denoted by 1-step-Intrpl. The reason 
for including 1-step-Intrpl is as follows: We choose a set 
of interpolation points and tangential directions; then 



using the interpolatory projection in Theorem 7, wc ob- 
tain a reduced model. Then, we use the same interpola- 
tion points and directions as an initialization technique 
for IRKA-PH. Then, we are able to illustrate starting 
from the same interpolation data, how IRKA-PH cor- 
rects the interpolation points and directions, and how 
much it improves the H2 and Hoo behavior through out 
this iterative process. 

Using all three methods, we reduce the order from r = 2 
to r = 20 with increments of two. For IRKA-PH, initial 
interpolation points are chosen as logarithmically spaced 
points between 10~^ and 10~^; and the corresponding 
directions are the dominant right singular vectors of the 
transfer function at each interpolation point (a small 
2x2 singular value problem). These same points and 
directions are also used for 1-step-Intrpl. The resulting 
relative H2 and ?^oo error norms for each order r are 
illustrated in Figure 3. Several important conclusions 
are in order: First of all, with respect to the both 7^2 
and "Hoc norm, IRKA-PH significantly outperforms the 
other methods. As the figure illustrates, for 1-step-Intrpl, 
the performance hardly improves as r increases unlike 
IRKA-PH for which both the 'H2 and "Hoo errors de- 
cay consistently. This shows the superiority of IRKA-PH 
over 1-step-Intrpl. The initial interpolation point and di- 
rection selection does not yield a satisfactory interpola- 
tory reduced-order model; however instead of searching 
some other interpolation data in an ad hoc way, the pro- 
posed method automatically corrects the interpolation 
points through out the iteration and yields a significantly 
smaller error norm. To see this effect more clearly, we 
plot the initial interpolation point selection, denoted by 
s^'^'^, and the final/ converged interpolation points, de- 
noted by s^*^"'*'^, in Figure 4 together with the mirror 
images of the original system poles, denoted by — Aj. As 
this figure illustrates, even starting with this logarith- 
mically placed points on the real line, IRKA-PH itera- 
tively corrects the points so that upon convergence they 
automatically align themselves around the mirror im- 
ages of the original systems poles. They don't necessar- 
ily exactly match the mirror images of the original poles 
but rather they align themselves in a way that overall 
they have a good representation of the original spectrum. 
This is similar to the eigenvalue computations where the 
Ritz values provide better information about the overall 
spectrum than just a subset of the exact eigenvalues. We 
further note that the full-order poles are computed only 
to obtain this figure and are not needed by IRKA-PH . 

The second observation from Figure 3 is that in com- 
parison to Eff. BaL, IRKA-PH achieves the better per- 
formance with less computational effort; the main cost 
is sparse linear solves and no dense matrix operations 
are needed unlike the balancing-based approaches where 
Lyapunov equations need to be solved. We also note that 
even though the proposed method is 'H2-based, it pro- 
duces a satisfactory Hoo performance as well. This is not 
surprising as [14] discusses IRKA usually leads to high 
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Relative H error as r varies Evolution of the relative error during IRKA-PH for r=20 




Relative error as r varies Evolution of the relative H_ error during IRKA-PH for r=20 




k: iteration index 



Fig. 3. Relative 'H2 and T-iao norms for Mass-Spring-Damper 
System 
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Fig. 4. Initial and converged interpolation points 
fidelity Hoo performance as well. 

Before we conclude this example, we further investigate 
the r = 20 case. We show, in Fig 5, how the ^2 and 
Hoo errors evolve during IRKA-PH for r = 20, the 
largest reduced-order. The figure reveals convergence 
within seven or eight steps. The large initial relative er- 
rors are reduced drastically already after two steps of 
the iteration. This illustrates the effectiveness of the pro- 
posed method. 

Next, we investigate the effect of different initializations 
on the performance of IRKA-PH. For brevity, we only 
illustrate the r = 20 case. We make 5 different initializa- 
tions and denote by Sq the set of initial interpolation 

points corresponding to the j*^ selection. Sq^^ will be 
the same as what was used earlier, i.e., 20 points loga- 
rithmically spaced between 10"'^ and 10^^. For Sq^\ we 



Fig. 5. Relative 'H2 and Hoo norms for Mass-Spring-Damper 
System 

choose 20 points logarithmically spaced between —10^^ 
and —10^^. Note that this is a poor selection since these 
initial interpolation points lie in the left-half plane in 
proximity of the poles. We consciously make this (bad) 

choice to see the effect on convergence. For S^^\ we 
choose complex points in the right half-plane with real 
parts that arc logarithmically spaced between 10~^ and 
10"' and imaginary parts that are logarithmically spaced 
between lO^'^ and 10"^. The set is chosen to be closed 
under conjugation. These points are arbitrarily selected 
and have nothing to do with the spectrum of A. For 
Sq^\ we make the situation even worse than for Sq^\ 
We choose 20 poles of the original system, G(s), and per- 
turb them by 0.1% to obtain our starting points. This 
is a very bad selection in two respects: 1) The interpo- 
lation points lie in the left-half plane; 2) They are ex- 
tremely close to system poles and make the linear sys- 
tem (siE— A)vi = Bbi very poorly conditioned. Finally, 

for Sq we choose, once again, r = 20 original systems 
poles, but this time reflect them across the imaginary 
axis to obtain initialization points. Once interpolation 
points are chosen, associated directions are taken to be 
the dominant right singular vectors of the transfer func- 
tion at each interpolation point as we did previously. 

Figure 6 shows the evolution of the relative 7^2 and Hoo 
errors during IRKA-PH for these 5 different selections. 
In all cases, IRKA-PH converges to the same reduced- 
model in almost the same number of steps. As expected, 

5o^'>and5^^> are the worst initializations, starting with 
relative errors bigger than 1 in the Hoo norm. How- 
ever, the algorithm successfully corrects these points and 
drives them towards high-fidelity interpolation points 

and tangent directions. Even though Sq seems to be 
the best initialization — it starts with the lowest initial 
error — it converges to the same reduced-model as the 
iteration that started with Sq^\ and also in the same 
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number of steps. Sq achieves this without the need for 
original poles. This numerical evidence is in support of 
the robustness of IRKA-PH. The algorithm successfully 
corrects the bad initializations. Indeed, IRKA-PH is ex- 
pected to be more robust and to converge faster than 
even the original IRK A, which already depicts fast con- 
vergence behavior. The reasons are twofold. Here, un- 
like in IRK A, regardless of the initialization, every in- 
termediate reduced-model is stable; hence the interpola- 
tions points never lie in the left-half plane. This already 
smoothens the convergence behavior. Moreover, unlike 
IRKA which uses an oblique projector, IRKA-PH is the- 
oretically equivalent to using an orthogonal projector 
(with respect to an inner product weighted by Q)) This 
in turn supports the robustness of the algorithm. 

For all these reasons, we expect IRKA-PH to be a 
numerically robust and rapidly converging iteration. 
Throughout our numerical experiments using a wide 
variety of different initializations, we have never expe- 
rienced a convergence failure of IRKA-PH. IRKA-PH 
always converged to the same reduced-model regardless 
of initialization; even a search for a counterexample 
spanning thousands of trials failed to produce even a 
single case of either convergence failure or convergence 
to a different model in the r — 20 case. Indeed this 
was true throughout the range 8 < r < 20; IRKA-PH 
converged to the same reduced-model for many differ- 
ent initializations for r = 8 : 2 : 20. Only for r = 2, 
r = 4, and r — 6 and then only after several trials, 
were we able to make IRKA-PH converge to a different 
reduced-model (only one different model) than what 
we had before. The performance of this new model was 
slightly better than what we had before; hence it would 
further support the success of IRKA-PH. However, as 
our simple initialization resulted in the same accuracy 
for every r = 8 : 2 : 20 and did only slightly worse for 
r = 2,4, 6, there is no need to update our earlier results. 
These numerical experiments once more strongly points 
toward the robustness of IRKA-PH. 

We present further comparison for the r — 20 case, 
by comparing time domain simulations for the full and 
reduced-models. Recall that G(s) has 2 inputs and 2 
outputs. To make the time-domain illustrations simpler, 
we only compare the outputs of the subsystem relating 
the first-input (ui) to the first-output (yi). The results 
for the other 3 subsystems depict the same behavior. 
As the input, we choose a decaying sinusoid Ui(t) = 
g-o.05t gjj-^ (^K^j-y j-^j^ ^j^g simulations for T = 50 sec- 
onds as it takes long time for the response to decay due 
to the poles close to the imaginary axis. The results are 
depicted in Figure 7. In Figure 7-(a), we plot the simula- 
tion results for the whole time interval. To give a better 
illustration, in Figure 7-(b), we zoom into the time in- 
terval [0, 10] seconds which illustrates that 1-Step-Intrpl 
leads to the largest deviation. In Figure 7-(c), we give 
the absolute value of the error between the true sys- 
tem output and the reduced ones. As it is clear from 
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Fig. 6. Relative 7^2 and Hao norms for Different Initializa- 
tions 



Full and reduced-models time domain simulations 
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Fig. 7. Time domain simulations for r = 20 

this figure and as expected from the earlier analysis as 
shown in Figure 3, IRKA yields the smallest deviation. 
The maximum absolute values of the output errors, i.e. 
max|yi(t) — yi,red(OI where yi,red(i) denotes the first- 
output due to the reduced-models are 
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6.2 MIMO port-Hamiltonian Ladder Network 

As a second MIMO port-Hamiltonian system, we con- 
sider an n-dimensional Ladder Network as shown in Fig. 



Relative error as r varies 
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Fig. 8. MIMO Ladder Network 



We take the current / on the left side and the voltage U 
on the right side of the Ladder Network as the inputs. 
The port-Hamiltonian outputs are the voltage over the 
first capacitor Uc-^ and the current through the last in- 
ductor II ,„ . The state variables are as follows: xi is the 
charge qi of Ci , X2 is the flux (pi of Li , is the charge q2 
of C2, Xi is the flux (j)2 of L2, etc. The directions chosen 
for the internal currents of the Network are shown by 
plus- and minus-signs in Fig 8. A minimal realization of 
this port-Hamiltonian Ladder Network for order n = A 
is given in Example 1. 

Adding another LC pair to the network, which would 
correspond to an increase of the dimension of the 
model by two, will modify the system matrices as 
follows: The subdiagonal of the matrix A will con- 



tain additionally i„/2-i' 



^ith the plus-sign in 



the (n/2 + l,n/2) position. The superdiagonal of A 
will contain — C~/„,L~,„ with the minus-sign in the 

n/2' n/2 o 

(rt/2, n/2 -I- 1) position. Furthermore, the main diagonal 

-R„/2 



of A will have 



2) position, 

-Rn/2+-Rn/2 + l 
' 1/2 



— in the (n — 2, n 

zero in the (n — l,n — 1) position, and 

in the (n, n) position. B and C matrices will be of the 
similar structure as in (15) and (16). The matrix B will 
have ones in the (1, 1) and (n, 2) positions and zeros in 
the other positions. The only nonzero entries of C are 
Yj- in the (1, 1) position and -j^ — in the (2,ri) position. 

1 71 / 2 

Here, we consider the 100-dimensional full order port- 
Hamiltonian network with Ci = 0.1 and Li = 0.1 and 
i?i = 3 for i = 1, . . . , 50 and R51 = 1. For this model, we 
compare LRKA-PH with both Eff. Bal. and regular bal- 
anced truncation (denoted by Reg. Bal). We would like 
to emphasize that Reg. Bal. does not preserve the port- 
Hamiltonian structure. We choose to include Reg. Bal. 
in our comparison to better illustrate the effectiveness 
of our proposed method by showing that IRKA-FH can 
perform as good as or sometimes even better than Reg. 
Bal. which is known to yield high-fidelity 'Hco and '7^2 
performance (though it is not constrained to preserve 
port-Hamiltonian structure). 




4 5 6 7 

Reiative H error as r varies 




Fig. 9. Evolution of the relative 'H2 and T-i^o norms 

We reduce the order from r = 1 to r = 10 in increments 
of one. The resulting relative H2 and T-Lco error norms 
are depicted in Figure 9. The ■H2-based nature of our 
proposed method is clear. IRKA-FH outperforms Eff. 
Bal. for every r. Indeed, what is interesting to observe 
is that IRKA-FH is better than even Reg. Bal. for each 
r = 1,...,5. Even though for r = 6,7,8 Reg. Bal. is 
better, for the last two r values, IRKA-FH is able to 
match Reg. Bal. Note that our proposed method achieves 
this performance while still preserving structure without 
any need for solving Lyapunov equations; the only cost is 
sparse linear solves. As expected, Reg. Bal. is the best in 
terms of T-Lao performance. Indeed, it is tailored towards 
H-oo error reduction and is not constrained to preserve 
structure. However, the 'H2-based IRKA-FH performs 
as well as Eff. Bal. in terms of "Hoc error norm. This once 
more shows that, similar to IRKA, IRKA-FH provides 
high-fidelity not only in the ^2 norm but also in the Hoc 
norm. 

Similar to the previous example, we illustrate the con- 
vergence behavior of IRKA-FH, in this case for r = 10, 
in Figure 10. In this case, the convergence is even faster 
than the previous example with the algorithm converg- 
ing after five to six steps. For this model, we have ini- 
tialized the interpolation points for IRKA-FH arbitrar- 
ily as logarithmically spaced points between 10~^ and 
10""^. As before, the initial tangent directions are chosen 
as the leading right singular vector of G(s) evaluated at 
the chosen interpolation points. 

As we did in the previous example, we experiment with 
different initialization strategies for IRKA-FH. For 
brevity, we only illustrate the r = 10 case. We consider 
5 different initializations, Sq for j = 1, . . . , 5 for the 

interpolation points. Sq^^ is what used earlier. For Sq^\ 
we choose r = 10 original systems poles and perturb 
them by 0.1% as the starting points. This is a very poor 
choice since the interpolation points are very close to 
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Evolution of the relative error during IRKA-PH for r=10 
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Fig. 10. Evolution of the relative "^2 and Hoo norms 



Evolution of the relative error during IRKA-PH for different initializatios 
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Fig. 11. Relative norms for Different Initializations 
131- 

the system poles. For Sq ,we choose 10 points logarith- 
mically spaced between — 10~^ and — 10~^. Once again, 
these interpolation points are also in the left-half plane; 
a strategy one would usually avoid. Sq^^ correspond 
to choosing 10 points logarithmically spaced between 

10-8 and 10-4. And finally, for S^^\ we choose some 
arbitrary complex numbers where the real parts are 
distributed between 10"^ and 10"^ and the imaginary 
parts are distributed between 10"^ and 10° together 
with their complex conjugate pairs. Once the starting 
points are chosen, the corresponding directions are the 
dominant right singular vectors of the transfer function 
at each interpolation point. Figure 11 shows the evolu- 
tion of the relative 7i2 error during IRKA-PH for these 
5 different selections. As before, in all five cases, IRKA- 
PH converges to the same reduced-model in almost 
the same number of steps illustrating the robustness of 
IRKA to to different initializations. 



6.3 Mass-Spring-Damper system with n — 20000 

In the previous two examples, to have a thorough anal- 
ysis of the models with as many system norm computa- 
tions as possible, we have chosen a very modest system 
size of 100. In this example, to illustrate that we can 
effectively apply our method in large-scale settings, we 
modify the Mass-Spring-Damper model to have a full 
model of order n = 20000. Then, using IRKA-PH, we 
reduce the order to r = 20, r — 30 r — AO and r — 50. 
In each case, the initial interpolation points are chosen 
as logarithmically spaced points between 10"'^ and 10"^ 
with the corresponding tangential directions chosen as 
the leading singular vectors of the transfer function at 
these points. Before we present the results, we note that 
even at this scale, the method took less than one minute 
to converge with a rather straightforward implementa- 
tion in Matlab. We have not tried to optimize the per- 
formance; we have simply used the Matlab sparse linear 
solves as is. The algorithm is expected to perform faster 
with the appropriate optimization of the code. 

The sigma plots, i.e. ||G(«cj)||2 vs. oj, of the full-order 
model G(s) and three of the four reduced models are 
plotted in Figure 12. We omitted the 40**^ order approx- 
imation to simplify the figure as the r = 40 approxima- 
tion was visually indistinguishable from the r — 50 one. 
Except the r = 20 case, all the reduced models provide 
a high quality approximation to the full-order model of 
order n = 20000. To illustrate the approximation qual- 
ity further, we depict the sigma plots of the error models 
in Figure 13. As r increases, the quality of the approxi- 
mation improves consistently; for r = 50, we obtain an 
approximatp~\ relative Hoo error of 7.90 x 10"^. Hence, 
with the proposed algorithm, we are able to reduce a 
port-Hamiltonian system of order n — 20000 in a numer- 
ically effective structure-preserving way using interpola- 
tion. Moreover, with the ■H2-inspired interpolation point 
and tangential direction selection, the reduced-model of 
order r = 50 is accurate to a relative error of 7.90 x 10-"*. 

The phase plots of the original model G(s) and the three 
of the four reduced-models are depicted in Figure 14. 
Once again we have left the 40"^ order approximation 
out due to the same reason as before. To make the figure 
simpler, we only depict the the phase plot of the subsys- 
tems relating the first input (ui(t)) to the first output 
(yi(<)). Similar to the sigma plot. Figure 14 illustrates 
that the plots for the full order model G(s) and for the 
reduced model G40 and G50 are virtually indistinguish- 
able. G30 shows some deviations around the low frequen- 
cies (as in the sigma plot). And as expected, G20 has the 



^ We use the term approsimaie since both ||G||koo and HG- 
GsoHwoo are computed by 500 frequency sampling points 
over the imaginary axis as opposed to an exact Hoc. norm 
computation. However, since the error plot is smooth, we 
expect this error number to be accurate enough. 
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^ Sigma plot for G{s) and the reduced-order models 
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Fig. 12. Sigma plots of G(s) and the reduced models 



Sigma plot for the error models 
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Fig. 13. Error sigma plots 

largest deviation. The phase plots for the other 3 sub- 
systems have the similar pattern and hence are omitted 
for the brevity of the paper. 

We conclude this example by illustrating time domain 
simulations for the same subsystem as above (from the 
first input Ui{t) to the first output yi{t)). We plot the 
results only for G(s), G2o(s) and 650(5) to make the 
illustrations more readable. First we use a decaying ex- 
ponential Ui — e~° °^* sin(5t) and run all three models 
for T = 50 seconds. The simulations are run using the 
Matlab ode45 solver (taking advantage of the sparsity 
of the matrices in the full-order simulations). The results 
are depicted in Figure 15-(a). Both reduced-models fol- 
low the full-order output very accurately. The maximum 
value of the absolute errors in the output responses are 
1.32 X 10-3 for G2o(s) and 6.15 x lO^^ for G5o(s). Ac- 
curate approximations are obtained in a fraction of the 
time, indeed. While the full-order simulation took 4.62 
seconds, the reduced-order simulations took 0.066 sec- 



Fig. 14. Phase plots of G(s) and the reduced models 

onds for G2o(s) and 0.068 seconds for G5o(s); more than 
98% reduction in simulation time. Of course, these gains 
will be significantly magnified when the full-order model 
needs to be simulated over and over again for different 
input selections. The second input we try is a square 
wave that oscillates between 1 and —1 with a period of 
0.27r seconds for the time interval — 20 second. Fig- 
ure 15-(b) shows the corresponding outputs. Once more, 
both reduced models depict a high-fidelity match. For 
this second input, the max value of the absolute errors 
in the output responses are 7.38 x 10^"^ for G2o(s) and 
5.85 X IQ-'^ for G5o(s). Once again, the simulation of the 
reduced-order models took a fraction of the time that 
took the full-order simulation. The simulation for G(s) 
took 9.15 seconds. On the other hand, it took 0.44 sec- 
ond for G2o(s) and 0.46 second for G5o(s); more than 
95% reduction in simulation time. As the original system 
dimension increases even further, these gains in simu- 
lation times would be one of the biggest advantages of 
model reduction. 
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8 Conclusions 

We have developed a framework for reducing multi- 
input /multi-output port-Hamiltonian systems via tan- 
gential rational interpolation. By choosing the projec- 
tion subspaces appropriately, we obtain reduced-order 
models that are not only rational tangential interpolants 
of the original system but they also retain the original 
port-Hamiltonian structure. Thus, they are guaranteed 
to be passive. An 'H2-inspired algorithm is introduced 
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Full and reduced-models time domain simulations 
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Fig. 15. Phase plots of G(s) and the reduced models 

for choosing the interpolation points and tangent direc- 
tions for high-fidelity approximations. Several numerical 
examples show the effectiveness of the proposed method. 
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